Time evolution of the reaction front in a subdiffusive system 
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Using the quasistatic approximation we show that in a subdiffusion-reaction system with arbitrary 
non-zero values of subdiffusion coefficients, the reaction front xj{t) evolves in time according to the 
formula Xf(t) = Kt a ^ 2 , with a being the subdiffusion parameter and K which is controlled by the 
subdiffusion coefficients. The parameter K is determined by the equation derived in this paper. To 
check correctness of our analysis, we compare analytical functions derived in this paper with the 
results obtained numerically for the subdiffusion-reaction equations. 
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I. INTRODUCTION 

The diffusion-reaction system with two initially sepa- 
rated diffusing particles of spices A and B reacting ac 



medium, as in porous media or gels [14|, [l5(. The 
subdiffusion is characterized by a time dependence of 
the mean square displacement of transported particle 
(Ax 2 ) = 2D a t a /T (1 + a), where D a is the subdiffusion 



cording to the formula m'A + n'B POnertl has been coefficient measured in the units m / s a and a is the sub 
intensively studied during past years [J, 0, 0> EL IE IE 
0, 0, 0, EE El El El- As the diffusion-reaction equa- 



tions describing the system are nonlinear, it is difficult to 
solve them and their general solutions remain unknown 
(except of very special cases). Thus, to simplify the cal- 
culations one usually uses various approximations, such 
as the quasistationary approximation [H, 0, 0] , the scal- 
ing method [l], 0, IE IE 0] > or the perturbation one [E EE] • 
Using these methods, there were derived characteristic 
functions of the system which include: the time evolu- 
tion of the reaction front Xf(t), the width of the reaction 
zone Wvdt) or the width of the depletion zone Wr, ep (t) 
QUE EL IE IE] which all appear to be the power functions of 
time f(t) = KV . The results were confirmed by numeri- 
cal calculations and simulations [1,0,0- However, as the 
methods of extracting the power functions are not based 
on analytical solutions of subdiffusion-reaction equations 
(not even on their approximately forms) the proportion- 
ality coefficient K is unknown. The coefficient carries 
dynamical information about the system e.g. how the 
diffusion coefficient influences the process. As far as we 
know, there were only a few attempts to determine of K 
by means of the quasistationary approximation 0, 0] . 

The situation is more complicated in the case of the 
subfiffusion system since the equations describing the 
system contain a derivative of fractional order. Subd- 
iffusion occurs in systems where mobility of particles is 
significantly hindered due to internal structure of the 
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diffusion parameter which obeys < a < 1. For a = 1 
one deals with the normal diffusion. Since no explicit 
solutions of the nonlinear (sub) diffusion-reaction equa- 
tions are known, one commonly considers a simplified 
system, for example a one in which diffusion coefficients 
of both reactants are assumed to be equal to each other. 
There is assumed that the some characteristic functions 
are the same in the system with any simplified assump- 
tions [IEE3- 

The problem is to choose a method to study the 
subdiffusion-reaction equations. The scaling method 
dose not allow one to determine K unless special condi- 
tions are taken into account. The perturbation method 
is of small efficiency because the first order correction is 
often insufficient, while the higher order corrections are 
hard to obtain even in the case of normal diffusion. So, 
there are a few problems solved using this method 9, 10]. 
An alternative method is the quasistationary one. In the 
case of normal diffusion-reaction system it is based on the 
assumption that process proceeds so slowly that changes 
of concentration of transported substance are small in 
some regions [E 0] ■ Since the subdiffusion process is sig- 
nificantly slower than the normal diffusion one, we expect 
that the quasistationary approximation is also applicable 
to the subdiffusive case. Thus, we adopt the method in 
this study. The scaling method and the quasistationary 
approximation one are often treated as equivalent to each 
other. We note however that the equivalence holds only 
in the long time limit [l7| . At shorter times applicability 
of the quasistationary method does not imply applica- 
bility of scaling one and vice versa (this problem will be 
discussed in fla|). 

In this paper we find that the time evolution of the 
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reaction front is given by the formula Xf (t) — Kt a / 2 for 
a system with arbitrary non-zero values of the subdif- 
fusion coefficients. The coefficient K fulfills the special 
equation derived in this paper. Our analytical results 
are confirmed by the numerical solutions of subdiffusion- 
reaction equations. 



II. THE SYSTEM 

A real system is usually three-dimensional, but we as- 
sume that it is homogeneous in the plane perpendicu- 
lar to the x axis. Therefore, we involve only one space 
variable x into considerations. The subdiffusion-reaction 
equations are 



d_ 

dt 



Ci{x,t)=D, 



d 1 



d 2 



dt 1 - dx 2 



C l (x,t) - diR a (x,t), (1) 



where i = A,B, denotes the concentration of the dif- 
fusing particles of species i, Di - the subdiffusion coef- 
ficient, cLa — to, ds — n, the parameters m and n oc- 
cur in the reaction term (Eq. ([3]) below); the Riemann- 
Liouvillc fractional time derivative is defined for the case 
of < a < 1 as 



d a f(t) 
dt a 



1 



r(i - a) dt 



f(r) 



dr- 

o (* " r) a 



Throughout this paper we assume that both of the reac- 
tants are mobile Da, Dg > 0. Let us note that the choice 
of the reaction term is not obvious [13, 0, H3, HH, [H, 
[23j . The reaction term, which we involve into consid- 
erations and which was used to study the subdiffusion- 
reaction system in [l6|, [ItJ , is 



R a (x,t) 



d 1 



R{x,t), 



(2) 



where the term R(x, t) within the mean field approxima- 
tion reads 



R(x,t) = kC%(x,t)C%(x,t), 



(3) 



k is the reaction rate and the parameters m and n are 
determined experimentally. 

We assume that the particles of reactants A and B 
are initially separated from each other. Thus, the initial 
conditions are 



C A {x,0) 
C b (x,0) 



Coa, x < 

0, x > 

0, x < 

Cos, x > 



(4) 
(5) 



It was observed 0, d, H, 0, [!, Q that when the process 
starts, there appear characteristic regions (see Fig. [T]): 
the depletion zone 'Dep', which is defined as a region 
where the concentrations are significantly smaller than 
the initial ones (Ca Co a and Cb <C Cob), the reaction 



region where the production of particles P is significant 
(R(x,i) > 0), and the diffusion region 'Dif, where the 
reaction term R(x, t) is close to zero and the particle 
transport appears to be almost subdiffusive (i.e. without 
chemical reactions). 




FIG. 1: Schematic view of the system under considerations; 
Xf(t) is the reaction front, 'Dep' and 'Dif' denote the deple- 
tion zone and the diffusion region, respectively. 

For the normal diffusion the widths of the depletion 
zone Wr> ep and the reaction region Wr grow as the power 
functions of time [H, 0, i, H H, H, 0| , W Dep ~ t 8 , with 
= 1/2, and W R ~ t", where (j, < 9. The value of 
parameter \x depends on the system under study. For the 
system where the reactants A and B are mobile, there is 
fi = 1/6 and for the system with a mobile reactant A and 
a static reactant B we have fi = (m— l)/2(m+ 1), where 
to is the parameter occurring in the reaction term R in 
Eq. © [| (see also 0,113). As was reported in 0, Il7j , 
Wr evolve in time according to the power functions also 
for the subdiffusive-reaction system with fi = a/6. 

An important characteristics of the system under con- 
sideration is the time evolution of the reaction front 
Xf(t). It is defined as a point where the reaction term 
R(x,t) reaches its maximum R(xf(t),t) = max or, as 
argued in for R <~ CaCb it is defined by the relation 
Ca (xf(t),t) = Cb(x f(t),t) and in more general situation 
by CA(xf(t),t)/m = C B (xf(t),t)/n 0. Unfortunately, 
these definitions are difficult to apply for the numerically 
obtained concentrations. In the following, we use the 
definition of the reaction front as 



*/(*) 



J xR(x, t)dx 
J R(x, t)dx 



(6) 



Although the relations defining the reaction front are not 
always equivalent to each other, all of them provide to 
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Xf lying inside the reaction region, and in the long time 
limit the definitions lead to the power functions of time. 
For the normal diffusion there is the dependence [J 0, 

naa 



(7) 



with 7 = 1/2. It was shown in [17| by means of the 
scaling method that the relation ((7|) with 7 = a/2 holds 
for the subdiffusive system where the subdiffusion coeffi- 
cients of reactants are equal to each other. 



III. QUASISTATIC APPROXIMATION 

The quasistatic approximation assumes that the con- 
centration profile is a slowly varying function of time in 
a given region. Thus, the time derivation is small and 
consequently, the r.h.s. of (sub)diffusion equation is 
also small in the region. It requires 



d 1 



d 2 



dt x - a dx 2 



C A: B{x,t) sa R a (x,t). 



(8) 



Since the reaction term is relatively large in the re- 
action zone, the quasistatic approximation holds in this 
zone under the condition 

D W^d^ CA ' Bix ' t] y> dl CAAx ' t] - (9) 

We note that the condition © is fulfilled when the con- 
centration profiles are given in the scaling form [l7| ■ 

In the diffusive region where R a (x,t) ss 0, the qua- 
sistatic approximation is applicable when the concentra- 
tion is a linear function of the r.h.s. ofEq. iQ} then 
vanishes. The regions outside the reaction zone, where 
the concentration linearly varies with x, determine the 
borders of the quasistatic region. The solution of subd- 
iffusion equation without chemical reactions works here. 

In the studies of the normal diffusion with reactions, 
one introduces the quasistatic approximation referring to 
the equilibration time tf @, S El ■ The parameter rp is 
of order of the average time which is needed for the sub- 
stance to spread over the interval of length Wr when 
the substance flows from outside of the interval. For the 
normal diffusion-reaction system this parameter was esti- 



mated from the relation (A2; 2 ) 



and t ~ tf one gets rp 
tern the relation (Ax 2 ) <~ t° 



t. Taking (Ax 2 ) 



For the subdiffusive sys- 
provides 



Tp - W, 



2/a 



(10) 



As for the normal diffusion case, let us assume that the 
relative change of the flux J fulfills the relation dJ/J = 
dt I tj which gives 



(tj)- 



rf(logJ) 
dt 



(11) 



The balance between the subdiffusion term and the reac- 
tion one is achieved when the equilibration time rp of the 
reaction region is negligibly small comparing to the time 
(tj) of relative change of the flux in the long time limit. 
So, the quasistatic approximation is applicable when 



0. 



(12) 



The quasistatic region is usually defined as a region where 
at least one of the conditions ([8]), ([9]) or (fT2]) is fulfilled. 
As far as we know, the equivalence of these definitions 
have not been proven yet. In our considerations we use 
the relation (jHJ) as the definition of quasistatic approxi- 
mation and we further show that the conditions (|12|) is 
fulfilled when Eq. (jHJ) is assumed. 



IV. TIME EVOLUTION OF Wr AND W Dep 

As in the normal diffusion reaction system, to find the 
widths of appropriate region we assume that the param- 
eters p — Db/Da and q — Cob /Co a are the irrelevant 
parameters of the system. This means that the results 
obtained for p = 1 and/or q = 1 are qualitatively equiv- 
alent to the one with p / 1 and/or q ^ 1, except of very 
few obvious cases (for example, when p = q = 1 the reac- 
tion front does not move). In this section we derive the 
time evolution of the widths of the reaction region Wr 
and the depletion zone Woep under condition p = q = 1 . 

At first we argue that the assumption \x < 6 is correct 
not only for the diffusive but for the subdiffusive systems 
as well. In [TtJ there was found that 9 = a/2 and p, = a/6 
by means of the simplified scaling method. We confirm 
the above relation, using the method already applied to 
the normal diffusion-reaction system (l3| |. 




FIG. 2: The symmetrical system with static reaction front 
xf(t) — 0. The continuous lines denote the concentrations C 
for the system under considerations, the dashed one - for the 
system with the fully absorbing wall Cabs located at x — 0. 

Let us assume that the initial concentrations and sub- 
diffusion coefficients of both reactants are equal to each 
other Cqa = Cob = Co and Da = Db = D. Due to 
the symmetry of the system, the reaction front will not 
change its position Xf(t) = 0. Proceeding as for the sys- 
tem with the normal diffusion [l3| , we assume to further 
simplify of the calculations that the concentrations of A 
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and B particles can be given as 

C A (x,t) = C Aahs (x,t) + 5C A (x,t), 

C B (x,t) = C Bahs (x,t) + 8C B (x,t), 

where C Aa bs and C Ba hs are the solutions of the pure 
subdiffusivc equation in the system with the perfectly 
absorbing wall located at x — and SC A and SC B 
are the corrections (see Fig. Symmetry of the sys- 
tem ensures that C A (x,t) = C B (—x,t), which provides 
SC A (x,t) = 5C B (—x,t). For a perfectly absorbing wall 
placed at x — 0, the concentration profiles vanish at the 
wall C Aa b s (0, t) = C Ba b s (0,t) = 0. After calculations, we 
obtain 



C Aabs (x,t) 



Co 



1 



10 



ttW 
- H ll 



2/a 



1 1 

2/a 



(13) 



for x < 0, and 

C Ba hs{x,t) 



1 




1 1 

2/a 



(14) 



for x > 0, where H denotes the Fox function, which can 
be expressed as the series [24| 



H 



1 1 

P q 



1 



q T- 



-IV 



f^ j!r(l - p/q - j/q) 



(15) 

Substituting C(x,t) = Ci(x.t), 8C(x,t) = 5Ci(x,t), 
C a bs(^, t) = Ciabsix, t), where i — A for x < and i = B 
for x > 0, to the subdiffusion-reaction equation and tak- 
ing into account that C a b s fulfills the subdiffusion equa- 
tion without chemical reactions, we get 



0_ 

di 



SC(x,t) 



d 1 



D^SC(x,t) 



dt 1 - 01 

-k [C ahs {x, t) + 6C(x, t)} SC{x, t)] 



The limits of the reaction region occur for x where 
SC is close to zero. In this region one can neglect 
the term (SC) 2 in the above equation. Moreover, in 
the long time limit we can approximate the Fox func- 
tions present in Eqs. 113[) and (fT4|) by the expression 
C ahs (x,t) = a\x\/t a / 2 , where a = C k/T(l - a/2)*/D. 
So, we get 



d_ 

di 



SC(x,t) 



d 1 



dt 1 



D JP § C(x,t)-^6C(x,t) 

OX t 2 



_ (16) 

As in the normal diffusion-reaction system [13j , we as- 
sume that 



0_ 

di 



8C(x,t) = 0. 



(17) 



From Eqs. ([16]). p7| and the relation [2 

dPt v r(z/ + i 



dtP r(v + 1 - (3) 



V > -1, 



we obtain 



D 



—&-5C(X,t) = -—; 
t 2 t a 



(18) 



(19) 



where A(x) is the arbitrary function of x only. In the long 
time limit the r.h.s. of Eq. (|19p can be neglected. Let 
us note that it is another justification to equal the r.h.s. 
of above equation to zero. To determine the function 
A[x) we observe that SC has a significant value in a finite 
region limited by the depending on time points —g(t) and 
g(t), which lie inside of the depletion zone (see Fig. 
Thus, we have 



SC(-g(t),t)^8C(g(t),t)^0 



and 



d 2 8C{x,t) 



dx 2 



d 2 SC(x,t) 



--git) 



dx 2 



--g(t) 



Since the left hand side of the Eq. (fl9|) is close to zero 
for |x| > g(t), the additional boundary conditions are 



A(-g(t)) = A(g(t)) = 0. 



(20) 



The above equations appear to be the boundary condi- 
tions for the function A, which cannot depend on time (in 
contrary to the boundary conditions (|20|1 ). Thus, there 
is the only solution A(x) = const = 0. 

Solving the equation ((19)) with the right side equal to 
zero, we find 



SC(x,t) = f(t)Ai (A-^), 



f 

where Ai denotes the Airy function, which can be ap- 
proximated by the following expression for large u 



Ai(u) 



1 



exp 



2u 3 / 2 



To obtain the function f(t) we assume that it is the power 
function of time f(t) ~ t x . Putting the function / to 



Eq. (fT6|) and using (fTTjl . we obtain A 
Eq. (H]) with (TT]) and ©, we get 

R(x,t) = ^P-SC(x,t). 

t 2 

Substituting Eq. ([21]) to Eq. dHJ) we obtain 

, 3 / 4 



Comparing 



(21) 



exp 



Ax 

t a/6 



3/2 



(22) 
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As the width of the reaction region is defined by the 
relation Q 



Wl(x,t) 



J(x - x f (t)) 2 R(x.t)dx 
J R{x, t)dx 



(23) 



it is easy to see that substituting (|22|) to (|23|) with Xf = 0, 
we obtain W R ~ i a/6 . 

Since the width of the depletion zone is defined by the 
conditions C t < C oi , i — A,B, from Eqs. (JT3J) and (H]) 
we get Wbcp ~ t a / 2 . Thus, the relation p < is fulfilled 
for the system where the subdiffusion coefficients of the 
reactants are equal to each other. We assume that this 
relation holds for the system with any non-zero values of 
the subdiffusion coefficients. 



V. CONCENTRATION PROFILE IN Dif 
REGION 

Since Wr <§C Wr> e p, the reaction region plays a role of 
a partially absorbing wall with respect to the depletion 
zone. We find the concentration profiles in the region 
outside the reaction one as a solution of the subdiffusion 
equation in the system with partially absorbing wall. 

We find here the solutions of the subdiffusion equation 
without chemical reactions (Eq. ([T]) with R a (x,t) = 0) 
for the system with partially absorbing wall. To calculate 
the concentration profiles, we use the integral formula 



C(x, t) = / G(x, t; x )C(xo, 0)dx o , 



(24) 



where G(x,t;xo) denotes the Green's function for the 
subdiffusion equation. From a macroscopic point of view, 
the Green's function is interpreted as a concentration 
profile of the N particles (divided by N) which are in- 
stantaneously produced and start from the position Xq 
at an initial moment t = 0. It is also interpreted as a 
probability density of finding a particle in a point x at 
time t under the condition that the particle is located in 
the position xq at the initial moment t = 0. 

There is a problem to set the boundary conditions 
at the partially absorbing wall. To obtain the Green's 
function one can use the method of images. The stan- 
dard method of images has been applied for the diffusive 
system with fully absorbing or fully reflecting wall [26j |. 
Then, one replaces the wall by a fictitious instantaneous 
point source of the particles (IPS) in such a manner that 
the concentration profile generated by all IPS behaves 
as in the system with the wall. In the system with the 
fully reflecting wall, the flux vanishes at the wall. In this 
case the Green's function can be obtained by replacing 
the wall by the auxiliary IPS of the same strength in the 
position symmetric to the initial point xq with respect to 
the wall 



where Go denotes the Green's function for homogeneous 
system. In the case of fully absorbing wall the concentra- 
tion vanishes at the wall. The Green's function is then a 
difference of IPS placed at xq and — xq, which gives 



G(x, t; x a ) = G (x, t; x ) - G (x, t; -x Q ). 



(26) 



Sometimes the boundary conditions are not given explic- 
itly by an equation, but they are postulated in a heuristic 
form. In such a case there is a possibility to use the gen- 
eralized method of images to find the Green's functions. 
Such a procedure was used to find the Green's functions 
for the system with partially permeable wall [27} where 
the Green's function was obtained from Eq. ([23)) by re- 
ducing the IPS located at — xq by the factor controlled 
by the permeability of the wall. 

For the system with partially absorbing wall we start 
with a physical condition, which can be stated as: if 
during a given time interval N particles reach the wall, 
the fraction p of them will be absorbed while 1 — p will go 
through. The parameter p is assumed to be a constant 
characterizing the wall. Such a situation appears when 
the partially absorbing wall is simulated by another IPS 
of the strength reduced by a factor p. So, the Green's 
functions are as follows 

GADif(x,t; xo) = G a{x, t;x ) - pAG A(x,t; -x Q ), 

(27) 

and 



G_BDif(a;, t;xo) = G b(x, t;x ) - p B G B(x,t; -x ), 



where 

G oi (x,t;xo) 



(28) 



a\x — xq\ 



H ll 



\X - Xq\ 




for i = A,B. Using the integral formula (|2"4"|) and initial 
conditions ^ and ([5]), we find (for details of the calcu- 
lations see the Appendix A) 



where 



and 



CAVit(x,t) = C a T)A 

a 

xH™ ( f^= X ~ 
11 \ \VDaT c 



Va = C A0 (1+pa)/2, 



Cb mf (x,t) = Cob rj B 

a 

2/a 




(30) 



(31) 



,VDb1^ 



1 1 

2/a 



(32) 



where 



G(x, t; x Q ) = G (x, t; x ) + G (x, t; -x a ), 



(25) 



Vb = C BQ {1 + p B )/2. 



(33) 
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Let us note that when Co a — Cqb — Cq and Da = Dg, 
we obtain pa = 1 and pg = 1 from Eqs. (fTB")) and (JH]). 
The subdiffusive fluxes are given by the formula 



Ji(x,t) 



dd(x,t) 



l dt 1 ~ a dx 
Using Eqs. |[30]) and (|32|). we obtain 



(34) 



4-i 



.VDa~^ 



2/a 



1 1 

-1 + 2/a 2/a 



(35) 



■TsDifOM) = \ZDbT]b ( — 



2/a 



1 



-1 + 2/a 2/a 



(36) 



In the following we use the shorter notation for the fluxes 
13 and 



J ATM 



t 1 -^ 2 * \^Da^J ' 
x 



JsDif - - t l- a /2 « 



where 



Q(*) 



2 ^ fe! r( a (i - k)/2) 



(-*)' 



(37) 
(38) 

(39) 



VI. TIME EVOLUTION OF THE REACTION 
FRONT 

In this section we derive the time evolution of the reac- 
tion front within the quasistatonary approximation. The 
derivation is based on three assumptions, which are ex- 
pected to hold in the long time limit. 

(1) We assume that the characteristic functions evolve in 
time according to the formulas 

W R ~ t a l\ 



(40) 



Bp 



«/2 



x s {t)~t«l\ (42) 

The relations were derived in [l7| by means of the 
scaling method for the system where the subdiffusion 
coefficients of both reactants are equal to each other. 
The relation (|4T)|) was also found in [8J| by means of the 
Monte Carlo simulations. Le us note that in Sec. IIVI 
we have shown that the relations (|40p and (|4"Tj) are 
fulfilled for the system where p — q = 1 . The relation 
will be also confirmed a posteriori in this section. 



(2) The region, where the quasistatic approximation 
works, extends beyond the reaction zone. Therefore 
there is the region defined by the relation 



W R (t) <\X- Xf(t)\ < WbepW, 



(43) 



where the quasistatic approximation region overlaps 
with the diffusion one. 

(3) In the diffusion region the concentrations are given 
by Eqs. (|37)l) - (f3"3")) with the parameters pa and pg, 
which can be larger than unity. 

Starting with the above assumptions, we show at first the 
following 

(a) The concentration profiles ([3H]) and (f3"2"| extended to 
the reaction region vanish at the points which are 
identified with the point x z defined in Fig. ([T]) and 
by Eq. (|47p . In the long time limit the point x z is 
localized so close to Xf that x z can be replaced by 

Xf. 

(Jo) The fluxes J a and Jb flowing into the reaction region 
from the left and from the right side, respectively, 
are assumed to be balanced in such a way that m 
particles A and n particles B flow into the reaction 
region in the time unit. 

After showing that the conditions (jlj) and Jb| hold, we 
use Eqs. (j30"f and (|32| to derive a relation describing the 
time evolution of the reaction front. 

As mentioned earlier, we are guided by the procedure 
already used for the normal diffusion-reaction systems 
. The quasistatic state can be defined by the following 
equations 



8 



l-a 



Of 1 



and 



8 



l-Q 



which combined provide 

Ql-a Q2 



d 2 

Da-^-jCa{x, t) — mR(x, t) 



d 2 

D b -q^Cb(x, t) - nR(x, t) 



0. 



= 0, 



dt x - a dx 2 



*(z,t) = 0, 



(41) where 



*(a;,t) = —D A C A {x,t) - -D B C B (x,t). 
m n 

Using the formula (fT8|) . we find 

t) = E{x)r a + F(t)x + G(t). 



(44) 



(45) 



Applying the operator 



|- to Eq. ([45]). we obtain 



d_ 

Of 



j^F(t) = -J B (t) - -J A (t). 



(46) 
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The function \& changes its sign in the reaction zone from 
positive where Cb — to negative where Ca — 0. Thus, 
there is the point x z (t) which lies inside the reaction zone, 
where the function ^ is equal to zero. Therefore, 

${x z (t),t)=0. (47) 

Since Xf(t) also lies inside the reaction zone there is 

\x z (t)-x f (t)\<nt a/6 , (48) 

where f2 is a positive constant. After simple calculations, 
we get 

T/ . E(x) — E(x z (t)) ... 

*0M) = ta +F(t)(x-x z (tj). (49) 

Let us now consider the region where the region of dif- 
fusion approximation overlaps with the one of the qua- 
sistatic approximation for x < xt(t). The region occurs 
for such x that the condition 

- W Bcp (x,t) <^x-x f {t) < -W R (x,t), (50) 

is fulfilled. Here Ca ~ Cabu, Cb ~ 0, J a ~ ^ADif, and 
J B ~ 0. So, we get from Eq. dHJ) 



V(x,t) = —D A C A r>u(x,t), 



m 



and from Eq. ([46]) 



F{ t) = — J Am t). 



eft 



(51) 



(52) 



Let us note that '5 is given by the function of the variable 
x/t a / 2 only (see Eq. ([30]) ). Therefore we deduce that 



^(a;) = ax 2 



(53) 
(54) 



Combining Eqs. ([30]), d49|), ([51]) and ®-([56]) we obtain 
e(t)-x f (ty 



D 



2 

Q 



Coa VaH u 



VDaI^ 





1 1 11 




2/a J 



a {x S {t)~e{t)f~xl(t) Xf (t)-e(t)-x z (t) 

a t a t a l 2 ' 1 ' 

Since in the longtime limit (xf(t) — x z (t) — e{t))/t a / 2 — > 
(see Eqs. (j48j and (57])), from Eq. |58]) we get 



1 1 

2/a 



= 0. (59) 



Similar considerations performed in the region 
W R (x, i) < x - x f (t) < W Dcp (x, t) 

provide 

1 



and 



y(x,t) = D B C B mi{x,t), 

n 



Ql-a i 

— F(t) = - JflDifW, 



dt 



which gives 
2 



c -. H io(( -fit) V " 
Cob - -VbH u ^ 



1 1 

2/a 



(60) 



0. (61) 



From Eqs. ([52]) and (|60[) we obtain 



1 

— Jadh 
m 



—Jbdh, 

n 



(62) 



and from Eqs. ([37]), ([35]) and §3§ we get 

1 



and 



(63) 

Combining Eqs. ([59]) . ([BT]) . ([63]) and using the identity 



G(t) = c, 



where a, 6, c are constants. 
We denote 



(55) 



Hl° ( z 2 ' a 
we have 



1 1 

2/a 



= 2 H n I * 



1 a/2 
1 



(64) 



x f (t)-x = e(t). (56) 

It is obvious that 

r>it Q/6 < e(t) < 2 t Q/2 , 

where Qi and ^2 are positive constants. When i — * oo 
the inequality provides W 6 /e(i) — > and 



i(t)/t a / 2 -> 0. 



(57) 



\y/DA~r) m^/D^CoB \VDb~^ 
1 a/2 



(65) 



where $(z) = ifjf z 



/Q(z). It is clear that 



1 

there is only one point Xf which for given t fulfills the 
definition of reaction front. The solution of the Eq. (IB51) 
is 



Xf (t) = Kt a / 2 , 



(66) 
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where coefficient K is the solution of the following equa- 
tion 



-K 



K 



(67) 



Thus, the time evolution of the reaction front is the power 
function with the exponent depending on the subdiffusion 
parameter a only; the subdiffusion coefficients Da and 
Db controll the parameter K . Eqs. (f6"6")l and (|6~T|) are 
the main result of our paper. 

The procedure developed in this paper is a extension 
of the one already used for the normal diffusion case Q. 
Repeating our consideration for a = 1 we obtain the 
results identical with those from ||. Our formula (|B"6"|) 
with K given by Eq. (|6"7| is a generalization of Eq. (21) 
in E3. 



VII. NUMERICAL SOLUTIONS 

To verify correctness of our procedure, we compare the 
analytical functions which are derived in the previous sec- 
tions with numerical solutions. We show that there exists 
the quasistatic approximation zone where, as required by 
Eqs. and the function is parabolic with 

respect to x. We also show that there exists the region 
of overlap of the diffusion zone and the quasistatic one; 
in this region Cadh or CsDif are the linear functions of 
x. 



A. Numerical procedure 

As we show in the Appendix B, assuming that the 
functions Ca and Cb and their second derivatives with 
respect to the space variable are limited, Eq. (JTJ) is equiv- 
alent to 

' da -Ci{x,t) = D~a(x,t)-d t R(x,t), (68) 



dt a 



dx 2 



where i = A, B, cLa = Tn, ds — n, with the Caputo 
fractional time derivative, which is defined for < a < 1 
as 29] 



d a f(t) 
dt a 



1 



r i 







Throughout this paper we denote the Riemann-Liouville 
fractional derivative without any additional index as 
d a f(t)/dt a , others kinds of the fractional derivatives are 
labeled by indexes C for the Caputo fractional derivative 
and GL for the Griinwald-Letnikov one. 

In the papers [3(1 HH there were presented the proce- 
dures of the numerical solving of the subdiffusion equa- 
tion without chemical reaction, when one can use the 
equation with Riemann-Liouville as well as Caputo frac- 
tional time derivative. The situation is different in the 



case of the subdiffusion-reaction equations. In the nu- 
merical calculations the fractional derivative is replaced 
by series. In the case of Eqs. ([T]) and {2} with Riemann- 
Liouville fractional derivative we have relatively compli- 
cated expression under the derivative, whereas in Eq. (|68|) 
the fractional derivative acts only on concentrations. It 
caused that the numerical procedure based on Eq. (|68|) 
is more convenient to use, at least in our opinion. 

To numerically solve the normal diffusion equation one 
usually substitutes the time derivative by the backward 



difference 



9 fit) 
dt 



f(t)-f(t-At) 
At 



In the presented proce- 



dure we proceed in a similar way. We use the Griinwald- 
Letnikov fractional derivative which is defined as a limit 
of a fractional-order backward difference 



dt° 



where a > 0, [z] means the integer part of z and 



(69) 



r(q + l) 
r J ~ r\T(a - r + 1) 

q(q - l)(a - 2) • . . . • [a - (r - 1)] 
1 • 2 • 3 ■ . . . • r ' 

When the function f(t) of positive argument has 
continuous derivatives of the first order, the Riemann- 
Liouville fractional derivative is equivalent to the 
Griinwald-Letnikov one for any parameter a (0 < a < 1) 
[23| . So, we have 



d a f(t) _ GL d a f(t) 



dt a 



dt a 



(70) 



The relation between Riemann-Liouville and Caputo 
derivatives is more complicated and reads as 



d a f{t) c d a f(t) 



dt c 



dt a 



+ #l_a(*)/(0), (71) 



where 



q+ w ' o t<o 



(72) 



From Eqs. ([69 |) - ([72|) we can express the Caputo fractional 
derivative in terms of the fractional-order backward dif- 
ference 

= ^)-^(-ir( a r )f(t-rAt) 



df 



t a T(l - a) 



/(0). 



(73) 



The standard way to approximate of the fractional 
derivative, which is useful for numerical calculations, is 
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to omit the limit in Eq. 
series to the finite one 



73|) and to change the infinite 



d a f(t) 
dt a 



(At)"* £)(-!) 



f(t - rAt) 



1 



t a T(l ~ a) 



/(0), 



(74) 



where the memory length L is a natural number of arbi- 
trary chosen value less then (or equal) [t/At]. 

Subdiffusion is a process with the memory. According 
to the short memory principle, the fractional derivative 
is approximated by the fractional derivative with mov- 
ing lower limit t — L, where L is the 'memory length' 
equals to a certain amount of time steps [29]. However, 



in the paper [32| there was shown that the numerical so- 
lutions of subdiffusion equation with the boundary con- 
ditions {J} and ([5|) are in agreement with the analytical 
one only when the memory length is closed to actual ac- 
count of time steps contrary to 'short memory principle'. 
So, in numerical calculations we take the memory length 
L equals to the actual number of time steps t s . 

Substituting Eq. ((74|) to Eq. ([68]). using the following 
approximation of the second order derivative 

d 2 f(x) f(x + Ax)-2f(x)+f(x-Ax) 



dx 2 



(Ax) 



we obtain 



d(x,t) 



E(-i 



r a(a- l)(a- 2) 



1-2-3 



t a T(l - a) 



+D 



(At) a 



[d(x + Ax,t- At) - 2Ci(x, t - At) + d{x -Ax,t- At)} 



(Ax) 2 

dik[At) a C%(x, t - M)C%{x, t - At), 



r 



(75) 



for i — A, _B, dA — rn and ds 



B. Numerical results 

Here we compare the analytical results with the nu- 
merical ones. In all figures there are presented functions 
calculated for the system where a = 0.5, Da = 0.025, 
D B = 0.0125, C A = 2, Cob = l,k = l,m = n= 1. For 
numerical calculations we take Ax = 0.2 and At = 0.05 
(all quantities are given in the arbitrary units). Addi- 
tionally, in Figs. ([5]) and ([6]) we plot the borders of the 
reaction zone (xf — Wr/2,x/ + Wr/2) calculated for the 
time 5000. The position of the reaction front was calcu- 
lated from the discrete version of Eq. ^ 



J2 t XtR(xj,t) 

J2iR{xi,t) ■ 



(76) 



and equals to 0.71 for t — 5000. The width of the reaction 
region calculated from discrete version of Eq. (f2"3")) 



wl(t) = 



Eiixj -x f {t)) 2 R{xj,t) 



equals to 0.38 for t = 5000. Thus the reaction region 
occupies the interval (0.52; 0.90). 
From Eq. (|76)) we find that 



This relation is very close to the relation 
calculated from Eq. (|67|) which reads 



x f (t) = 0.0825t 



0.25 



with K 



(78) 



x f (t) = 0.0838t 



0.251 



(77) 



In Figs. [3] and |4] there are presented the concentration 
profiles Ca and Cb obtained numerically according to 
the formula (|75|) and the functions given by Eqs. (|3"0|) 
and (|3"2")1 with pa = 0.40 and ps = 3.64, respectively. 
We observe a quite good agreement of the analytical and 
numerical functions in the diffusion region. 

In Fig. O we present the function ^(x, t) calculated 
numerically and its parabolic approximation ^S(x,t) = 
0.297(x/t Q / 2 ) 2 - 0.168(x/t Q / 2 ) + 0.015. We note that * 
is satisfactorily approximated by the parabolic function 
of x. The region where ^ is parabolic determines the 
quasistatic approximation region. 

In Fig. [5] we present the numerical solutions of 
the subdiffusion-reaction equations and their linear ap- 
proximations calculated from the formulas CU(x,t) ~ 
-0.816a; + 0.616 and C B (x,t) m 0.620x - 0.490, respec- 
tively. As seen, the linear approximation of Ca and Cb is 
satisfactory outside the reaction region. This statement 
confirms correctness the quasistationary approximation 
in the region enclosing the reaction region. 

We conclude this section by saying that our numerical 
results support the postulates of the quasistatic approx- 
imation. 
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FIG. 3: The symbols represent the numerical solutions of 
subdiffusion-reaction equation, the continuous lines are as- 
signed to theoretical functions Cadh for the times given in 
the legend. 




FIG. 4: The symbols represent the numerical solutions of 
subdiffusion-reaction equation, the continuous lines are as- 
signed to theoretical functions Cboh for the times given in 
the legend. 

VIII. FINAL REMARKS 

Using the quasistationary approximation and utilizing 
the solution of the subdiffusion-reaction equations in the 
diffusive region, we show that the time evolution of the 
reaction front for the subdiffusion-reaction system is a 
power function (|66[) with the exponent a/2 and the coef- 
ficient K is controlled by the subdiffusion coefficients of 
the system. The function Xf ~ t a ^ 2 can be obtained by 
means of the scaling method. However, it is very hard 
within this method to find an explicit expression of the 
parameter K for the case of Da ^ Db- 




X 



FIG. 5: The function 9 (symbols) obtained numerically for 
the times given in the legend and their parabolic approxima- 
tions inside the quasistatic approximation region (continuous 
line); the vertical lines represent the borders of the reaction 
zone calculated for t = 5000. 



2,0-. 




1,6- 
1,4- 
1,2 - 



0,8- 

0,6- 




-5 -4 -3 -2 -1 1 2 3 4 5 
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FIG. 6: The concentration profiles Ca and Cb obtained nu- 
merically (squares) calculated for time 5000 and their linear 
approximations (dashed lines), the vertical lines represent the 
borders of the reaction zone. 



We note that in this paper we consider the process of 
subdiffusion controlls chemical reactions. It means that 
the reactions which proceed relatively fast when com- 
pared to the characteristic time of meeting of the par- 
ticles of A and B [l[ . Under such assumption the qua- 
sistatic approximation works and the time evolution of 
the reaction front does not depend on the detailed form 
of the reaction term (expect of dependence of the pa- 
rameters m and n). This happens because the form of 
R does not change the relation Wn ~ i Q / 6 . Thus, the 
width of the reaction zone appears to be relatively small 
in comparison with the width of the quasistatic approxi- 
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mation region. The time evolution of Xf is determined by 
the dynamics of transport of the particles to the reaction 
zone and it depends on the parameters m, n, Da, Db, 
Cqa and Cos only. This statement is particularly im- 
portant for the subdiffusion-reaction systems where the 
reaction term is not uniquely defined (as the fractional 
derivative can be involved into this term in a few ways 

EE E3, El HI EH, IH 

As far as we know, the time evolution of the reaction 
front has not been measured experimentally in a subdif- 
fusive system with two mobile reactants. For this reason 
we can compare the functions ([50)) and (f3"2"]> with exper- 
imental data obtained for a subdiffusive system without 
chemical reactions. Our theoretical and the experimen- 
tal functions presented in [3| are qualitatively similar to 
each other if we take the units commonly used in real 
systems where x is given in 10 _2 m, t in sec, Da and 
Db are of the order 10 _8 m 2 / s a . Since K is controlled 
by the subdiffusion coefficients of reactants, the method 
presented in this paper can be used for extracting the 
subdiffusion parameter from experimental data. The nu- 
merical calculations show that if we take the subdiffusion 
coefficients of the order maintained above and we assume 
that Cqa/Cob is of the order of 1, we obtain K ~ 10~ 2 
m/s a/2 from Eq. ([67)1 . 

The quasistatic approximation in a normal diffusion 
system applies to a region where the equilibrium time 
Tp of the reaction region is negligibly small comparing to 
the characteristic time of change of the flux tj in the long 
time limit 0, 0, EH ■ Let us note that this fact is fulfilled 
in the subdiffusive-recation system. Since Wr ~ i Q / 6 , 
we have tf ~ t 1 / 3 form (| 10[) . Taking the definition 



(jTTj) . which for the subdiffusion flux J ~ 1/i 1 Q / 2 gives 
tj - 1/t (see Eqs. (37]) and (38])), we get t f /tj 
for any value of the subdiffusive parameter a. So, the 
assumptions adopted in our paper agree with the quasis- 
tationary condition Q12]). 

The function <3> is approximated by parabolic function 
in the region where the quasistatic approximation region 
overlaps with the diffusion one. However, we expect that 
there are departures from this approximation in a re- 
gion located within the reaction zone where the reaction 
term is significantly different from zero. It is because 
of the concentrations C a and Cb have different scaling 
properties in that region. We expect that the width of 
that region is so narrow, as compared with the width the 
quasistatic approximation one, that the departure form 
parabolic approximation is hard to observe on the plots 
presented in our paper. We note that the possibility of 
occurring this departure does not influence our main re- 
sults. 
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APPENDIX A 

In this Appendix we present some details of the pro- 
cedure of solving the subdiffusion equation. The calcula- 
tions with Ricmann-Liouville fractional time derivative 
are relatively simple in terms of the Laplace transform 
(LT) L[f(t)\ = f(s) = J °° dtf(t)e- st . The LT of Green's 
function for a homogeneous subdiffusive system without 
chemical reaction (|29[) reads [IH, H3| 



G ,i(x,s;x ) 



1 



2jD~s 1 - a / 2 



-\x-x \^s"/Di^ (Ai) 



i = A, B. The LT commutes with the integration (whith 
respect to the variable x), so from Eq. (fM| we get 



C(x,s) = / G(x, s; xo)C(xo,0)dxo. 



(A2) 



Putting (jATj) to the LT of Eqs. (27]) and (28]). and next 
to Eq. (|A2|) . we obtain 



(A3) 



(A4) 



C A (x,s) = °™ - 3de-(— Q/2 )/v^. 
s s 

with r/A = Cqa(1 + Pa)/2, and 



s s 



withes = Cos(l+ps)/2. The inverse Laplace transform 
L^ 1 gives (here a > and > 0) [H[ 



L 
and 



- 1 (>e— ") 



1 



L- l (s v e- aa " 



^ oo 



ttU) 



,l//3 



1 1 

l + v 1 



(A5) 



(A6) 

Using the relation (|A5[) to calculate the inverse LT of 
Eqs. O and (H]) we get (3D]) and (3JJ). Let us note 
that comparing the right hand sides of Eqs. (|A5[) and 
(|A6[) . after simple calculations we get the useful relation 

The LT of the subdiffusive flux (3JD reads 
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Ji(x, s) = -DiS 



i- a dCi{x,s) 
dx 



(A7) 



Putting Eqs. (A3| and (|A4|) to ([AT]) and next to Eq. ([A5]) . 
we obtain ([33]) and (31)]) . 
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APPENDIX B 

Here we show that Eqs. ([lj and ([68]) are equivalent 
to each other when the concentration C and its second 
space derivative are limited. The Laplace transforms of 
fractional derivatives are as follows [29| (here < a < 1) 



L 



d a f(t) 
dt a 



l f(t) 



dt° 



t=o 



d a f(t) 
dt a 



= s a f(s)-s a - 1 f(0). 



The Laplace transform of Eq. ((T|) is 



sC(x,s) ~ C(x,0) = 
d 2 C(x,s) 



D- 



dt- a 



D 



dx 2 

d 2 C(x ,t) 
dx 



R(x, s) 
, -R(x,t) 



(Bl) 



whereas the Laplace transform of Eq. (|68p reads as 

'■d{ 

d.l 



s a C{x,s) - s Q - 1 C(a;,0) = D * c k' s) - R(x,s), which 



gives 

sC(x,s) — C(x, Q) = s 1 - 01 



D 



d 2 C{x,s) 
dx 2 



i?(x, s) 



(B2) 



We assume that the function C and its second space- 
variable derivative are limited. So, there is a positive 
number M which fulfils the relation \<d(x, t)\ < M, where 

Q(x, t) = D d2 ^f t] -R{x, t), for any (x, t). From the def- 
inition of Riemann-Liouville derivative of negative order 



dt- 



we obtain 



:/(*) = 



1 



r (") Jo 



dr(t-T) a -V(r), 



-G(x,t) 



df- 



and from Eq. (|B3j) we get 



<M f dr(t - t) 
Jo 



a 



dt 



^eOM)| t=0 = o, 



(B3) 



which causes that the Laplace transforms (|T31|) and (|B2|) 
are equal to each other. Thus, the equations |T]) and 
are equivalent to each other for the limited function 9. 
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